####Replication for Bove and Gavrilova (2017) -----
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
#This file replicates Tables 1-4 from B+G, as well as Figures 1 and 2.
#To run, you need the militarization.dta and police.dta dataframe.
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
library(tidyverse)
library(readstata13)
library(xtable)
library(plm)
library(lfe)
library(texreg)
library(stargazer)
library(haven)
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
#CLEAN DATA AND CREATE VARIABLES OF INTEREST ------
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
####load first data frame and calculate variables
df <- read.dta13("raw-data/militarization.dta")

###all of this is to calculate the instrument - see page 8 of B+G
df <- df %>%
  mutate(ind = ifelse(weapons_quantity!=0 | 
                        vehicles_quantity!=0 |
                        gears_quantity!=0 |
                        other_quantity!=0,1,0),
         ind = ifelse((year<2006), NA, ind),
         ind = if_else(is.na(ind), 0, ind)) 

df <- df %>% 
  group_by(fips) %>% 
  summarise(Aid1 = sum(ind))%>%
  left_join(df,.) %>%
  arrange(fips,year) 

BG_milit <- df %>%
  group_by(fips) %>%
  mutate(burden_iv = Aid1/7*dplyr::lag(burden,2),
         milex_iv = Aid1/7*log(dplyr::lag(milex,2)),
         milex_iv_3yr = Aid1/3*log(dplyr::lag(milex,2)),
         casualties_iv = Aid1/7*log(dplyr::lag(casualties,2))) %>%
  select(-statetimefe_2,-statetimefe_3,-statetimefe_7, ##drop vars dropped by stata
         -statetimefe_133,-statetimefe_139,-statetimefe_156,
         -statetimefe_159,-statetimefe_161,-statetimefe_166,
         -statetimefe_174,-statetimefe_176,-statetimefe_181,
         -statetimefe_182,-statetimefe_184,-statetimefe_185) 

save(BG_milit, file = 'data/BG_milit.RData')

# these variables cause problems in Stata because they have missingness
temp <- BG_milit %>%
  select(-state,-county)
write.dta(temp, file='data/BG_milit.dta')


# create restricted data frame: 2010 to 2012
BG_milit_restricted <- BG_milit %>%
  filter(year>2009)
  # select(-statetimefe_2,-statetimefe_3,-statetimefe_7, ##drop vars dropped by stata
  #        -statetimefe_133,-statetimefe_139,-statetimefe_156,
  #        -statetimefe_159,-statetimefe_161,-statetimefe_166,
  #        -statetimefe_174,-statetimefe_176,-statetimefe_181,
  #        -statetimefe_182,-statetimefe_184,-statetimefe_185)

#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
#TABLE 1 ----------
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
# Defining the function
my.summary = function(x, na.rm=TRUE){
  result = c(Mean=mean(x, na.rm=na.rm),
             SD=sd(x, na.rm=na.rm),
             Median=median(x, na.rm=na.rm),
             Min=min(x, na.rm=na.rm),
             Max=max(x, na.rm=na.rm), 
             N=length(x))
}

# identifying numeric columns
ind = c(61:67, 39:44, 261,262, 68:69, 5:12, 14:15, 17, 52, 54:59)

##drop any observations that are NA (via calculation of lags) to match sample
##size of 17822
summarydf <- BG_milit %>%
  filter(!is.na(milex_iv) & !is.na(crime)) 

summarydf.res <- BG_milit_restricted %>%
  filter(!is.na(milex_iv) & !is.na(crime))  

# applying the function to numeric columns only
table1 = as.data.frame(sapply(summarydf[, ind], my.summary))  
table1 = t(table1)

table1 <- as.data.frame(table1)
table1$Variable <- c('Crime rate','Murder rate','Robbery rate',
                     'Assault rate','Burglary rate','Larceny rate',
                     'Vehicle theft rate','Murder arrest rate',
                     'Robbery arrest rate','Assault arrest rate',
                     'Burglary arrest rate','Larceny arrest rate',
                     'Vehicle arrest theft','Military exp. Iv','Military exp. Iv adj.',
                     'Total aid (value)','Total aid (quantity)',
                     'Weapons aid','Weapons quantity','Vehicles aid',
                     'Vehicles quantity','Gears aid','Gears quantity',
                     'Others aid','Others quantity','Percent poverty',
                     'Median income','Unemployment rate','Population',
                     'Share males','Share blacks','Share age 15–19',
                     'Share age 20–24','Share age 25–29','Share age 30–34')
row.names(table1) <- NULL
table1 <- table1 %>%
  select(Variable,Mean,SD,Median,Min,Max,N)
colnames(table1)[7] <- "Observations"

out.table1 <- xtable(table1,
                     caption = "Summary Statistics",
                     digits=c(0,0,1,1,1,1,1,0))
print(out.table1,
      booktabs = TRUE,
      include.rownames = FALSE,
      file = "tables/BGtable1.tex",
      caption.placement = "top",
      format.args = list(big.mark = ",", decimal.mark = "."))
# table 1 for restricted data --------------
# applying the function to numeric columns only
table1.res = as.data.frame(sapply(summarydf.res[, ind], my.summary))  
table1.res = t(table1.res)

table1.res <- as.data.frame(table1.res)
table1.res$Variable <- c('Crime rate','Murder rate','Robbery rate',
                         'Assault rate','Burglary rate','Larceny rate',
                         'Vehicle theft rate','Murder arrest rate',
                         'Robbery arrest rate','Assault arrest rate',
                         'Burglary arrest rate','Larceny arrest rate',
                         'Vehicle arrest theft','Military exp. Iv','Military exp. Iv adj.',
                         'Total aid (value)','Total aid (quantity)',
                         'Weapons aid','Weapons quantity','Vehicles aid',
                         'Vehicles quantity','Gears aid','Gears quantity',
                         'Others aid','Others quantity','Percent poverty',
                         'Median income','Unemployment rate','Population',
                         'Share males','Share blacks','Share age 15–19',
                         'Share age 20–24','Share age 25–29','Share age 30–34')
row.names(table1.res) <- NULL
table1.res <- table1.res %>%
  select(Variable,Mean,SD,Median,Min,Max,N)
colnames(table1.res)[7] <- "Observations"

out.table1.res <- xtable(table1.res,
                         caption = "Summary Statistics (2010-2012)",
                         digits=c(0,0,1,1,1,1,1,0))
print(out.table1.res,
      booktabs = TRUE,
      include.rownames = FALSE,
      file = "tables/BGtable1-restricted.tex",
      caption.placement = "top",
      format.args = list(big.mark = ",", decimal.mark = ".")) 
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
#TABLE 2 -----------
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
#first, get naive OLS estimates for COLUMN 1
fols = as.formula(paste('crime ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                        lpop +ltotal_cost+', 
                        paste(colnames(BG_milit)[54:59], collapse='+'), 
                        '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|0',
                        '|State'))

felm1<-felm(fols, data=BG_milit) 
summary(felm1)

#then COLUMNS 2 AND 3 (FIRST STAGE OF IV AND REGULAR RESULTS OF THE IV)
fa = as.formula(paste('crime ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                      lpop +', 
                      paste(colnames(BG_milit)[54:59], collapse='+'), 
                      '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|(ltotal_cost~milex_iv)',
                      '|State'))

felm2<-felm(fa, data=BG_milit)

summary(felm2)

#COLUMN 4
fhomicide = as.formula(paste('murder ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                             lpop +', 
                             paste(colnames(BG_milit)[54:59], collapse='+'), 
                             '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|(ltotal_cost~milex_iv)',
                             '|State'))


felm3<-felm(fhomicide, data=BG_milit)

summary(felm3)

#COLUMN 5
frobbery = as.formula(paste('robbery ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                            lpop +', 
                            paste(colnames(BG_milit)[54:59], collapse='+'), 
                            '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|(ltotal_cost~milex_iv)',
                            '|State'))

felm4<-felm(frobbery, data=BG_milit)

summary(felm4)

#COLUMN 6
fassault= as.formula(paste('assault ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                           lpop +', 
                           paste(colnames(BG_milit)[54:59], collapse='+'), 
                           '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|(ltotal_cost~milex_iv)',
                           '|State'))

felm5<-felm(fassault, data=BG_milit)

summary(felm5)

#COLUMN 7
fburglary= as.formula(paste('burglary ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                            lpop +', 
                            paste(colnames(BG_milit)[54:59], collapse='+'), 
                            '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|(ltotal_cost~milex_iv)',
                            '|State'))

felm6<-felm(fburglary, data=BG_milit)

summary(felm6)

#COLUMN 8
flarceny= as.formula(paste('larceny ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                           lpop +', 
                           paste(colnames(BG_milit)[54:59], collapse='+'), 
                           '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|(ltotal_cost~milex_iv)',
                           '|State'))

felm7<-felm(flarceny, data=BG_milit)

summary(felm7)

#COLUMN 9
fmvtheft= as.formula(paste('mvtheft ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                           lpop +', 
                           paste(colnames(BG_milit)[54:59], collapse='+'), 
                           '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|(ltotal_cost~milex_iv)',
                           '|State'))
felm8<-felm(fmvtheft, data=BG_milit)

summary(felm8)

# Table 2 Compile -----------------
#look at results - replicated successfully - FOR TABLE 2
name.map <- list(milex_iv = "Military exp. IV",
                 milex_iv.3yr = "Military exp. IV adj.",
                 ltotal_cost = "Lagged Total Aid",
                 "`ltotal_cost(fit)`" = "Lagged Total Aid")

model.list <- list(felm1,felm2$stage1,felm2,felm3,felm4,felm5,felm6,felm7,felm8)

texreg(model.list,
       file = "tables/BGtable2.tex",
       use.packages = FALSE,
       custom.coef.map = name.map,
       digits = 3,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = c("OLS","Firststage","Crime","Homicide","Robbery",
                              "Assault","Burglary", "Larceny","Vehicletheft"),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = F)

temp <- map(model.list,extract)
saveRDS(temp,"data/BG-table2-original.RDS")

# Table 2 with restricted data (2010 to 2012) -----------------

#first, get naive OLS estimates for COLUMN 1
fols.r = as.formula(paste('crime ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                          lpop +ltotal_cost+', 
                          paste(colnames(BG_milit_restricted)[54:59], collapse='+'), 
                          '| fips+',paste(colnames(BG_milit_restricted)[81:194], collapse='+'),'|0',
                          '|State'))

felm1.r<-felm(fols.r, data=BG_milit_restricted) 

#then COLUMNS 2 AND 3 (FIRST STAGE OF IV AND REGULAR RESULTS OF THE IV)
fa.r = as.formula(paste('crime ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                        lpop +', 
                        paste(colnames(BG_milit_restricted)[54:59], collapse='+'), 
                        '| fips+',paste(colnames(BG_milit_restricted)[81:194], collapse='+'),'|(ltotal_cost~milex_iv)',
                        '|State'))

felm2.r<-felm(fa.r, data=BG_milit_restricted)

#COLUMN 4
fhomicide.r = as.formula(paste('murder ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                               lpop +', 
                               paste(colnames(BG_milit_restricted)[54:59], collapse='+'), 
                               '| fips+',paste(colnames(BG_milit_restricted)[81:194], collapse='+'),'|(ltotal_cost~milex_iv)',
                               '|State'))


felm3.r<-felm(fhomicide.r, data=BG_milit_restricted)

#COLUMN 5
frobbery.r = as.formula(paste('robbery ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                              lpop +', 
                              paste(colnames(BG_milit_restricted)[54:59], collapse='+'), 
                              '| fips+',paste(colnames(BG_milit_restricted)[81:194], collapse='+'),'|(ltotal_cost~milex_iv)',
                              '|State'))

felm4.r <-felm(frobbery.r, data=BG_milit_restricted)

#COLUMN 6
fassault.r = as.formula(paste('assault ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                              lpop +', 
                              paste(colnames(BG_milit_restricted)[54:59], collapse='+'), 
                              '| fips+',paste(colnames(BG_milit_restricted)[81:194], collapse='+'),'|(ltotal_cost~milex_iv)',
                              '|State'))

felm5.r <-felm(fassault.r, data=BG_milit_restricted)

#COLUMN 7
fburglary.r = as.formula(paste('burglary ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                               lpop +', 
                               paste(colnames(BG_milit_restricted)[54:59], collapse='+'), 
                               '| fips+',paste(colnames(BG_milit_restricted)[81:194], collapse='+'),'|(ltotal_cost~milex_iv)',
                               '|State'))

felm6.r <-felm(fburglary.r, data=BG_milit_restricted)

#COLUMN 8
flarceny.r = as.formula(paste('larceny ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                              lpop +', 
                              paste(colnames(BG_milit_restricted)[54:59], collapse='+'), 
                              '| fips+',paste(colnames(BG_milit_restricted)[81:194], collapse='+'),'|(ltotal_cost~milex_iv)',
                              '|State'))

felm7.r <-felm(flarceny.r, data=BG_milit_restricted)

#COLUMN 9
fmvtheft.r = as.formula(paste('mvtheft ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                              lpop +', 
                              paste(colnames(BG_milit_restricted)[54:59], collapse='+'), 
                              '| fips+',paste(colnames(BG_milit_restricted)[81:194], collapse='+'),'|(ltotal_cost~milex_iv)',
                              '|State'))
felm8.r <-felm(fmvtheft.r , data=BG_milit_restricted) 

model.list.res <- list(felm1.r,felm2.r$stage1,felm2.r,felm3.r,felm4.r,felm5.r,felm6.r,felm7.r,felm8.r)

texreg(model.list.res,
       file = "tables/BGtable2-restricted.tex",
       use.packages = FALSE,
       custom.coef.map = name.map,
       digits = 3,
       caption.above = TRUE,
       caption = "The Effect of Military Aid on Crime",
       custom.model.names = c("OLS","Firststage","Crime","Homicide","Robbery",
                              "Assault","Burglary", "Larceny","Vehicletheft"),
       #dcolumn = TRUE,
       booktabs = TRUE,
       include.rsquared = F)

temp.res <- map(model.list.res,extract,include.rsquared = F,include.adjrs = F)
saveRDS(temp.res,"data/BG-table2-original-restricted.RDS") 

#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
#TABLE 3 -----
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
#COLUMN 1
fhomicide = as.formula(paste('amurder ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                             lpop +', 
                             paste(colnames(BG_milit)[54:59], collapse='+'), 
                             '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|(ltotal_cost~milex_iv)',
                             '|State'))

felm2<-felm(fhomicide, data=BG_milit)
summary(felm2)

#COLUMN 2
frobbery = as.formula(paste('arobbery ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                            lpop +', 
                            paste(colnames(BG_milit)[54:59], collapse='+'), 
                            '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|(ltotal_cost~milex_iv)',
                            '|State'))

felm3<-felm(frobbery, data=BG_milit)
summary(felm3)

#COLUMN 3
fassault= as.formula(paste('aassault ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                           lpop +', 
                           paste(colnames(BG_milit)[54:59], collapse='+'), 
                           '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|(ltotal_cost~milex_iv)',
                           '|State'))

felm4<-felm(fassault, data=BG_milit)
summary(felm4)

#COLUMN 4
fburglary= as.formula(paste('aburglary ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                            lpop +', 
                            paste(colnames(BG_milit)[54:59], collapse='+'), 
                            '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|(ltotal_cost~milex_iv)',
                            '|State'))

felm5<-felm(fburglary, data=BG_milit)
summary(felm5)

#COLUMN 5
flarceny= as.formula(paste('alarceny ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                           lpop +', 
                           paste(colnames(BG_milit)[54:59], collapse='+'), 
                           '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|(ltotal_cost~milex_iv)',
                           '|State'))

felm6<-felm(flarceny, data=BG_milit)
summary(felm6)

#COLUMN 6 
fmvtheft= as.formula(paste('amvtheft ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                           lpop +', 
                           paste(colnames(BG_milit)[54:59], collapse='+'), 
                           '| fips+',paste(colnames(BG_milit)[81:249], collapse='+'),'|(ltotal_cost~milex_iv)',
                           '|State'))

felm7<-felm(fmvtheft, data=BG_milit)
summary(felm7)

## Table 3 compile ---------------------
model.list <- list(felm2,felm3,felm4,felm5,felm6,felm7)
table3.out <- map(model.list,extract)
saveRDS(table3.out,file="data/BG-table3-original.RDS")
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
#PREPARATION FOR TABLE 4 ------
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
##FIRST, CLEAN DATA AND CONDUCT DATA MANIPULATION
#first, import police.dta and repeat the same data manipulation
dfpolice <- read.dta13("raw-data/police.dta") %>%
  select(-statetimefe_1, -statetimefe_9,-statetimefe_11,
         -statetimefe_13,-statetimefe_17,-statetimefe_19,
         -statetimefe_21,-statetimefe_23,-statetimefe_29,
         -statetimefe_31,-statetimefe_37,-statetimefe_41,
         -statetimefe_46,-statetimefe_48,-statetimefe_52,
         -statetimefe_166,-statetimefe_173,-statetimefe_192,
         -statetimefe_194,-statetimefe_197,-statetimefe_203,
         -statetimefe_212,-statetimefe_214,-statetimefe_220,
         -statetimefe_221,-statetimefe_223,-statetimefe_224)  #drop vars from dataset that Stata also drops

##select observations Stata is dropping to drop them from this dataset too
#(stata is dropping cases that only have one observation)
singletons <- dfpolice %>% 
  group_by(agency) %>%
  summarise(number = n())   

#now, merge dropped observations with original dataframe to delete ones with one 1 obs
df2 <- left_join(dfpolice,singletons) %>%
  filter(number>1) 

#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
#TABLE 4-------------------
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
#COLUMN 1
fixed_effects <- grep("statetimefe",names(df2))
fa = as.formula(paste('officerx1000pop ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                      lpop +', 
                      paste(colnames(df2)[44:49], collapse='+'), 
                      '| agency+',paste(colnames(df2)[fixed_effects], collapse='+'),'|(ltotal_cost~milex_iv)',
                      '|State')) 

felm1<-felm(fa, data=df2)
summary(felm1)

#COLUMN 2
femp = as.formula(paste('tot_empl ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                        lpop +', 
                        paste(colnames(df2)[44:49], collapse='+'), 
                        '| agency+',paste(colnames(df2)[fixed_effects], collapse='+'),'|(ltotal_cost~milex_iv)',
                        '|State')) 

felm2<-felm(femp, data=df2)
summary(felm2)

#COLUMN 3
fofftoemp = as.formula(paste('off_to_empl ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                             lpop +', 
                             paste(colnames(df2)[44:49], collapse='+'), 
                             '| agency+',paste(colnames(df2)[fixed_effects], collapse='+'),'|(ltotal_cost~milex_iv)',
                             '|State'))

felm3<-felm(fofftoemp, data=df2)
summary(felm3)

#COLUMN 4
ftotalcalls = as.formula(paste('total_calls ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                               lpop +', 
                               paste(colnames(df2)[44:49], collapse='+'), 
                               '| agency+',paste(colnames(df2)[fixed_effects], collapse='+'),'|(ltotal_cost~milex_iv)',
                               '|State'))

felm4<-felm(ftotalcalls, data=df2)
summary(felm4)

#COLUMN 5
finjuries = as.formula(paste('injuries ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                             lpop +', 
                             paste(colnames(df2)[44:49], collapse='+'), 
                             '| agency+',paste(colnames(df2)[fixed_effects], collapse='+'),'|(ltotal_cost~milex_iv)',
                             '|State'))

felm5<-felm(finjuries, data=df2)
summary(felm5) 

#COLUMN 6
fasslts = as.formula(paste('civil_disorder_asslts ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                           lpop +', 
                           paste(colnames(df2)[44:49], collapse='+'), 
                           '| agency+',paste(colnames(df2)[fixed_effects], collapse='+'),'|(ltotal_cost~milex_iv)',
                           '|State'))

felm6<-felm(fasslts, data=df2)
summary(felm6)

#COLUMN 7
foffenderskilled = as.formula(paste('offenders_killed ~ PovertyPercentAllAges + lMedian + Unempl_Rate + 
                                    lpop +', 
                                    paste(colnames(df2)[44:49], collapse='+'), 
                                    '| agency+',paste(colnames(df2)[fixed_effects], collapse='+'),'|(ltotal_cost~milex_iv)',
                                    '|State'))

felm7<-felm(foffenderskilled, data=df2)
summary(felm7)

#COLUMN 8
fcalls= as.formula(paste('logcom ~ PovertyPercentAllAges + lMedian + Unempl_Rate +
                         lpop +',
                         paste(colnames(df2)[44:49], collapse='+'),
                         '| fips+',paste(colnames(df2)[fixed_effects], collapse='+'),'|(ltotal_cost~milex_iv)',
                         '|State'))

felm8<-felm(fcalls, data=df2)
summary(felm8)

##check results FOR TABLE 4
model.list <- list(felm1,felm2,felm3,felm4,felm5,felm6,felm7,felm8)
table3.out <- map(model.list,extract)
saveRDS(table3.out,file="data/BG-table4-original.RDS") 
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
#FIGURE 1
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
tmp = aggregate.data.frame(df[, c(61, 72)], list(df$fips), function(x) mean(x, na.rm=TRUE))

plot(tmp$crime, tmp$total_raw_cap, main="Figure 1", 
     xlab="Crime Rate ", ylab="Total Aid per Capita ", xlim = c(0,10000),
     ylim = c(0,20))

rm(tmp)

#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
#FIGURE 2
#----------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------------------------------------#
df$Aid1 = ifelse((df$Aid1==6|df$Aid1==5|df$Aid1==4|df$Aid1==3), 7, 1)

tmp = aggregate.data.frame(df[, c(61, 70)], list(df$year, df$Aid1, df$milex), function(x) mean(x, na.rm=TRUE))
tmp = tmp[-c(1:4),]
tmp$Group.3 = tmp$Group.3/1000
tmp = tmp[order(tmp$Group.1),]
tmp7 = tmp[tmp$Group.2==7,]
tmp1 = tmp[tmp$Group.2==1,]

plot(tmp7$Group.1, tmp7$crime, type = 'l', main = 'Figure 2a', xlab = 'Year', ylab = 'Total Aid', 
     ylim = c(2000, 4000))
lines(tmp1$Group.1, tmp1$crime, col='blue')
par(new = TRUE)
plot(tmp7$Group.1, tmp7$Group.3, type = "l", col = 'red', axes = FALSE, bty = "n", xlab= '', ylab = "")
axis(side=4, at = pretty(range(tmp7$Group.3)))
legend("topleft",legend=c("High recipients","Low recipients", 'US military spending'),
       text.col=c("black","blue", 'red'),col=c("black", 'blue', "red"))

plot(tmp7$Group.1, tmp7$ltotal_cost, type = 'l', main = 'Figure 2b', xlab = 'Year', ylab = 'Total Aid', ylim = c(0, 8))
lines(tmp1$Group.1, tmp1$ltotal_cost, col='blue')
par(new = TRUE)
plot(tmp7$Group.1, tmp7$Group.3, type = "l", col = 'red', axes = FALSE, bty = "n", xlab= '', ylab = "")
axis(side=4, at = pretty(range(tmp7$Group.3)))
legend("topleft",legend=c("High recipients","Low recipients", 'US military spending'),
       text.col=c("black","blue", 'red'),col=c("black", 'blue', "red"))

rm(tmp)
rm(tmp1)
rm(tmp7)
